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Materials with Periodic Internal Structure:^ 
Computation Based on Homogenization and 
Comparison with Experiment 


S. Jansson, F.A. Leckie, E.T. Onat, and M.P. Ranaweera 


Abstract 


The combination of thermal and mechanical loading expected in practice means that 
constitutive equations of metal matrix composites must be developed which deal with time- 
independent and time-dependent irreversible deformation. Also, the internal state of 
composites is extremely complicated which underlines the need to formulate macroscopic 
constitutive equations with a limited number of state variables which represent the internal 
state at the micro level. One available method for calculating the macro properties of 
composites in terms of the distribution and properties of the constituent materials is the 
method of homogenization whose formulation is based on the periodicity of the 
substructure of the composite. 

In this study a homogenization procedure has been developed which lends itself to 
the use of the finite element procedure. The efficiency of these procedures, to determine 
the macroscopic properties of a composite system from its constituent properties, has been 
demonstrated utilizing an aluminum plate perforated by directionally oriented slits. The 
selection of this problem is based on the fact that, i) extensive experimental results 
exist, ii) the macroscopic response is highly anisotropic and iii) that the slits provide very 
high stress gradients (more severe than would normally be found in practice) which severly 
test the effectiveness of the computational procedures. Furthermore, both elastic and 
plastic properties have be investigated so that the application to practical systems with 
inelastic deformation should be able to proceed without difficulty. The effectiveness of the 
procedures have been rigorously checked against experimental results and with the 
predictions of approximate calculations. Using the computational results it is illustrated 
how macroscopic constitutive equations can be expressed in forms of the elastic and limit 
load behavior. 


t Work funded under NASA Grant NAG3-894. 
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INTRODUCTION 


The combination of thermal and mechanical loading expected in practice means that 
constitutive equations of metal matrix composites must be developed which deal with time- 
independent and time-dependent irreversible deformations. The so-called unified 
constitutive equations are likely to provide a good basis for the description of the matrix 
material but the combined effect of the matrix and reinforcing material remains to be 
determined. One available method for calculating the macro properties of composites in 
terms of the properties and distribution of the constituent materials is the method of 
homogenization whose formulation is based on periodicity of the substructure of the 
composite. The method can coincide with the classical method of representing the 
substructure by a unit cell representation when conditions of symmetry are valid. However 
the advantages of the method of homogenization are that its formulation allows it to be used 
when symmetry no longer applies, offers the possibility of determining the stress and strain 
fields at the microscopic level and provides a formal derivation of unit cell representation. 
By following this procedure failure criteria can be introduced into the calculations for the 
constituents and interfaces which are based on micro-mechanical models. Studies which 
establish the macro-mechanical properties of the composite from the properties of the 
constituents afford the opportunity of directly designing composite material properties. 
When studying component behavior however it is more convenient to use constitutive 
equations which describe the macroscopic properties of the material. Such constitutive 
equations can be developed from the results of mechanical tests. However it is also 
possible to use homogenization procedures to simulate the experimental program and 
provide the results from which constitutive equations can be formulated. 

In this study a convenient form of this homogenization procedure has been developed 
which lends itself to the use of the finite element procedure. The procedure is based on the 
assumption that stress fields vary slowly from one homologous point to another. In the 
absence of point loads and away from boundaries we expect this assumption to be valid. 
However, this expectation is not fulfilled in the vicinity of stress-free boundaries or in the 
vicinity of cracks where stress gradients are large. For nonlinear and inelastic problems 
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localization of deformation is a possibility. If localization takes place in a region of a few 
cells then homogenization may not be applicable because macro-strain gradients are large. 
However, the procedure may still be valid provided localization occurs over many cells 
(twenty say). 

To illustrate the application of the procedure a plate system has been selected for 
which extensive experimental results exist. The system consists of an aluminum plate 
perforated by directionally oriented slits. In addition to providing a basis of comparison for 
the predictions of the calculations, the stress concentration factors are high and provide a 
severe test on the effectiveness of the computational procedures. Finally the anisotropic 
character is pronounced and is difficult to represent at the macroscopic level. 

A further study given in a subsequent report is the prediction of the macroscopic 
properties of a metal matrix composite using the known properties of the matrix and the 
fiber. An experimental program has also been completed so that comparison with the 
computational results can be made. This study is particularly important since if the 
predictions of the calculations can be verified then the procedure can replace the difficult 
test program which is often limited by the shortage of material. 

2. THE HOMOGENIZATION PROCEDURE 

2.1 Formulation 

Consider a structure composed of a periodically inhomogeneous material (Fig. 2.1). 
The material is linear elastic and locally isotropic. For further simplicity assume that the 
inhomogeneity is planar so that the the Lam 6 constants of the material depend only on the 
transverse coordinates xj and X 2 as follows: 

(i (x + d e a ) = ji(x), x e R 2 , a =1,2 (2.1) 

where ej and e 2 are the first two unit vectors of the rectangular frame and x lies in their 
span (similar equation for the other Lame constant X). 

We shall be concerned here with small plane deformations and rotations which 
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in the present case can take place in the absence of the stress components Oj3 and ©23 and 
we shall regard the material as infinite and two dimensional. It is then useful to think of the 
structure as a two dimensional sub-set B of the material shown in Fig. 2.1 where the 
material is composed of the copies of a unit cell of square cross section which is of size d 
and contains, say, a hard fiber of circular cross section. 

The elastic structure B is subject to given surface tractions T on dBp and given 
displacements U on dBy. It is required to determine the resulting fields ui(x), efj(x), and 
CTj(x) of displacement, strain and stress respectively (Fig. 2.1). 

If the characteristic lengths of this elasticity problem (the diameter D of the body, the 
minimum radius of curvature of the concave parts of the boundary 0B of the body, the 
"wave" lengths 1 T and ly associated with surface tractions and displacements, etc.) are 
much larger than the cell size d, then one expects that the solution of the above problem of 
elasticity would exhibit certain properties. A careful statement of these properties will 
require a family of decreasing cell size so that statements can be made about the material 
properties as d/D tends to zero. Thus following the French school^) we will define the 
Lame constants of the material with cell size d as (i^\x) = jo.^)(x/d) where is a one- 
periodic continuous function. 

For a discussion of the anticipated properties of the elasticity problem of Fig. 2.1 it is 
desirable to consider a generic nine-cell sub-domain within the body (Fig. 2.2a). The 
deformed shape of this domain may be as shown in Fig. 2.2b. Due to inhomogeneity of 
the material, the stress fields within a cell will exhibit variations so that a(Pj), the stress at 
Pi. is likely to be quite different from o(Qi) and the deformed shape of the boundary of 
the cell will exhibit a waviness of size d. However, if one considers three homologous 
points Py P2 and P3 within the sub-domain (these points have the same coordinates in 
respective cell coordinate frames), then cr(Pi), cj(P 2) and o(P3) will differ from each other 
only a little. More precisely, say 

||a(P 1 )-a(P 2 )||/(llall)«l 
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In view of the above cited expectations, the deformed state of the nine-cell sub- 
domain can be a pproximated by the one shown in Fig. 2.2c which is macroscopically 
homogeneous. Thus, in this figure, the same stress will be obtained in homologous 
locations Pj' which now lie on a straight line and, therefore, the new locations A'B'C'C' of 
the cell-comers will define a parallelogram. 

Reflection will show that the displacement field Vj(x) associated with the 
macroscopically homogeneous field shown in Fig. 2.2c will be of the following form in the 
case of small plane deformations and rotations 

Vi(x) = Eij xj + Qij xj, ij = 1,2 (2.2) 

where the usual summation convention is used and Ejj and Qjj are the components of the 
constant infinitesimal tensors of strain and rotation that determine the shape and orientation 
of the parallelogram A'B'C'D' Here, Vi(x) are the components of a d- period ic and 
continuous displacement field, 

Vi(x + d ej) = Vi(x) ( 2 - 3 ) 

The strains created by the field (2.2) are 

eijW^Vjj+V+Eij (2.4) 


where comma denotes partial differential in the usual way. 

It will be noted that on account of the periodicity of Vj the average of eij over a cell is 
equal to Ejj: 

R=e« < 2 - 5 > 

where the brackets indicate the average, over a cell, of the bracketed quantity. 
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The stresses created by the above strains will satisfy the equations of equilibrium: 


_a_ 

8xj 


Cijkito 


2 


(V k ,i + V ljk ) 


= 0 


( 2 . 6 ) 


where Qjktfx) are the d-periodic moduli of elasticity 

Qjld 00 = |x(x)(5i k 5ji + 5ii 5j k ) + A.(x)5i jk i . 

It is known that, when combined with the periodicit y of continuity of Vj(x), Eqn. 
(2.6) constitutes a well-posed linear problem defined over the interior of a cell for the 
determination of Vj(x) as functions of the prescribed constant average strain components 

% 

In view of the linearity of the problem, and with the arbitrary but permissible choice 
that Vj(x) = 0 on the comers of cells, (2.6) has the unique solution 

Vi(x) = gij k (x)E jk (2-7) 

where the functions gjj k (x) follow uniquely from the elastic moduli C(x) of the material. 
Note that 


§ijk ®ikj 

and there are six such d-periodic functions. In this report the functions g ijk (x) are 
determined by finite element methods. 

Combining (2.7) with (2.4) and the Hooke's law, we find the macrostrains and 
stresses created by the constant average strain E k p 

Eij(x) = Ejj + l/2(gi k i, j(x) + gjlddOO) E kl> 
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where 


Gij(x) = Hjjki(x) Ey 


(2.8) 


HjjklOO “ [QjklOO + QjmnOO §mkl> n( x )l- 

From the expression (2.7) it can be deduced that are of the order d and therefore 
gmkl.n ^ of the order "one" [2,33- 

Note that Hjjki(x) are d-periodic. It follows from the equations of equilibrium and 
from the fact that Ey are arbitrary constants that 

Hijkl( x )>j = tQjkl( x ) + Qjmn( x ) 8mkl,n( x )]J = 0 (2.9) 

It will be observed that (2.9) could be thought as providing six equations for the 
determination of the function gmkl( x ) an d is equivalent to (2.6). 

By averaging, we obtain from (2.8) 

Zij=K ijld E kl (2.10) 

where 

Zy - (<Jy) c and K ijkl = (H ijkl (x)) c (2. 1 1) 

and brackets indicate, as before, the average of the bracketed quantity over a cell. 

We refer to Zjj as the average stress or the macro-stress. Kjjkl are the components of 
the constant effective moduli of the problem. It can be shown by using (2.9) that Kjjkl 
inherit the symmetry properties of Qjki exhibited before. Moreover it can also be shown 
that for the family of the materials considered here Kjjkl are independent of the cell size d! 

It is noted that finite element techniques can and have been applied (cf. [1]) to a single 
cell to obtain approximations to the functions gijk( x )> Hjjki(x) and, hence to the effective 


7 


moduli Kjjifj. 

After this initial step which is independent of the body and the loads it carries and 
produces the effective moduli and the important relationships (2.7) and (2.8), we are ready 
to discuss the remaining steps of the homogenization program which leads to an 
approximate solution of the problem of finding the fields of displacement and stress within 
the body B. 

2.2 Homogenized Body. Macro-Stresses and Strains 

The structure B is now thought of as composed of a homogeneous elastic material 
which possesses the effective moduli K obtained in Section 2.1. The distributions E(x) 
and E(x) of the macro- or average-stresses and strains are then found by solving a classical 
boundary value problem of elasto-statics over B with the boundary conditions of the 

original problem (Fig. (2.1): 

Ejj nj = Ti(x) on 3Bp 

Uj = uj(x) on 3By (2.12) 

where Uj(x) are the displacements that give rise to the distribution Ejj(x) of the macro- 
strains over the body. 

The distribution of the micro stresses and strains within the body can then be 
determined in the following way. 

In view of (2.7) and (2.2) the displacement field uj is offered as an approximation to 
the actual field uj(x) of the very first problem of Fig. 2.1): 


uj(x) = Uj(x) + gijk(x)E jk (x) 


(2.13) 


where Ui(x) and Ej k (x) are the fields of displacement and strain in the homogenized body 
and the functions gijkW were obtained in the first step of the homogenization program. 

The above displacement field is kinematically acceptable everywhere except on the 
boundary BBy where the second term on the r.h. of the above equation and therefore the 
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difference (u, - Uj) will, in general, not vanish (cf. Fig. (2.2b)). 

Let p“ denote the deficit in displacements on 8Bu: 

Pi = Uj(x) - Uj(x) = gjj k (x)Ej k (x) (2.14) 

We recall that gjjk is of the order of d and therefore we expect that the steps needed to 
be taken to remove this deficit will cause only local changes in the fields introduced above. 
Errors shall be discussed later. 

The strains created by the displacement field Uj(x) are then 

Gij( x ) = Ejj(x) + ^(gjjjjj + gjk],i)E k ](x) + l/2(gi k jE k jj + gjlcjEy j) (2.15) 

which give rise to the stress field 

djj(x) = C ijkl (x)e kl (x), xeB (2.16) 

In view of the developments which led to (2.8) and (2.11), the above equation 
becomes 


djj(x) — Hjj k ](x)E k j(x) + [Cjj k j(x)g kmn (x)E mn j(x)] (2. 17) 


where as before 

Hijkl( x ) = [Qjkl( x ) + QjmnOO gmkl,n( x )l- 

We note that the first term on the right hand side of (17) is of the "order" of the 
macro- stresses £jj(x) and coincides with the leading term of the asymptotic expansion used 
in the work of Duvunt [2] for the stress components. We also observe that the second term 
in (2.17) is of the order of d/D. 
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2.3. Errors in Homogenization 

It is a valid question to pose how close are the fields Uj(x) and Gjj(x) obtained by the 
homogenization program to those occurring in the heterogeneous medium. 

It was shown in [2] that 

Zjj(x) - HyuM E k i(x) — » 0 (2.18) 

as d/D — ■» 0. 

This important result does not, of course, imply that the above difference tends to 
zero everywhere. Indeed, as we shall see, on 3Bp the difference is not zero, in general, 
and remains so as d/D tends to zero. Moreover when d/D is finite but still small there may 
be other parts of the body where the above difference may be unacceptably large (cf. Fig. 
2. 2d). 

The starting point in the construction of the approximate fields was the choice of the 
displacements uj(x). It was demonstrated that this field fails to satisfy the b.c. on 3By and 
it is now studied whether the stress field d;j(x) derived from it satisfies the requirements of 
equilibrium. 

Consider first the equilibrium of the material elements adjacent to the boundary dBp 
of the body where surface tractions are prescribed. Equilibrium demands that the vector 

a i j(x)nj(x)-T i (x) 

should vanish on dBp where nj(x) are the components of the unit outward normal to the 
boundary. Since CTji(x) exhibits, in general, spatial variations over the cell due the 
inhomogeneity of the material and T[(x) is, by assumption (cf. Fig. 2.1), a function of long 
wave length the above difference will not, in general, be zero (cf. Fig. 2.2e). 

This leads us to introduce the vector value function p|(x) which measures the deficit 
in equilibrium on dBp: 
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Pi(x) = Oij(x)nj(x) - Tj(x) on 9 B f 


(2.19) 


Note that in view of (2.12) the above equation can be written as 

p[(x) = (a,j(x) - £ij(x))n j(x)) on 9B F (2.20) 

In the case shown in Fig. 2.3e max Ip 2 (x)l on the stress free boundary is independent 
of d and it is known from the numerical work reported later that it is of the order of the 
applied stress 5^ x . 

Similarly the deficit q(x) in equilibrium in the interior B of the body is defined by 
considering the divergence of the stress field djj(x): 

q(x) = ajjj(x) (2-21) 


and using (2.17) 

fi(x) = Hijjci(x),j E^i(x) + Hjjj^x) Ej c ](x),j 
+ (Qjkl( x ) 8kmn( x ) E mn,l( x ))]’j 

The first term on the r.h. of the above vanishes by (2.9) so that 

q(x) = H ijkl (x) Ej^xXj + [Cijki(x) gkmn( x ) E mn,l( x ))].j ( 2 - 22 ) 

We observe that each of the two terms on the right hand of (2.22) will be of the order 
of lldy(x)ll/D in the interior of the body B. However near a concave portion of the 
boundary 3B the gradient of Ej^x) can be large so that a more refined estimation of q(x) 
will be needed in such locations. 

It is important to note that the solution of the boundary value problem for the 
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inhomogeneous elastic body B with the boundary conditions 

Ui = -pFoo on 9B U 

T i = — Pi ( x ) on 3b F 

and subject to the body forces bj (2.23) 

bi = -ri(x) inB 

would provide the fields of displacement and stress which when added respectively to the 
fields uj(x) and ajj(x) would produce the exact fields uj(x) and <?ij(x). 

Clearly the' problem of estimating errors in homogenization is equivalent to that of 
estimating the solution of (2.23). This leads us to the study of the deficits listed on the 
right hand of (2.23). 

The specific example selected for illustration involves corrections for the deficit in 
equilibrium on 3B f . Consider first a reduced version of the elastostatics problem defined 
in (2.23). Let uj and bj vanish, but let Tj have the values given in (2.23). The question 
then arises about the displacements and stresses created in B by the surface tractions — pj(x) 
applied on 3Bp? The important result (2.18) of asymptotic analysis suggests that when 
d/D is small these correction fields must decay strongly from 3Bp towards the interior of 
the body. The reason for this strong decay is a version of the de Saint Venant principle 
when applied to periodically inhomogeneous elastic bodies. 

For this purpose consider the following integral 

Pi 

f(Pi.P2)= J (a u (x) - Z u (x))dx 2 , X 1 = P1 (2.24) 

P2 

where the segment of integration which is of length d and parallel to the x 2 -axis is assumed 
to be on the boundary of B. 
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Note that if the point (p^^) where on dBp and if Xj were constant on the boundary 
then (2.24) would measure the resultant of p[(x) over the side of a cell: 

p 2 +d 

f(Pl>P 2 )= J Pi(x)dx 2 , xj = P! (2.25) 

P2 

Now using (2.17) and (2.10), but neglecting terms of order d/D in (2.17), with the 
expansion 

E kl(Pl> x 2) " E kl(Pl-P2) + E kl,2(Pl>P2)( x 2 - P2> 

+V2!E kli22 (pi,P2)( x 2 " P2) 2 

(2.24) becomes 

P2+d 

f(Pl»P2) ss J (Hiiki(Pi.x 2 )-K ilk i)[E k i(p 1 ,p 2 ) + E k i >2 (pi,p 2 )(x 2 -p 2 ) 

P2 

+V2!E kl>22 (p 1 ,p 2 )(x 2 -p 2 ) 2 ]dx 2 (2.26) 

Using (2.9) and (2.1 1) and the d-periodicity of H(x) gives the result 

p 2 +d 

J (H 1 ikl(Pi. x 2 )-K i iki)dx 2 = 0 (2.27) 

P2 

so that 

p 2 +d 

f(Pl> P 2 ) * J ( H ilkl(Pl- x 2)- K ilkl)[ E kl,2(Pl'P2)( x 2 _ P2) 

P2 

+1/2 ! E k ] 22 (p 1 , p 2 )(x 2 - p 2 ) 2 ]dx 2 (2.28) 
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This is an important result on various counts. First of all it shows that 


f( Pl , P2 ) - rf-ujta.pilCXd 2 ) (2.29) 

where k is a constant which has the dimension of stress. 

In the case where the boundary has a section parallel to the X 2 *axis (2.29) and (2.25) 
show that the resultant of the deficit p[(x) over a length d of the boundary is tWQ orders o£ 
magnitude smaller than the local value of the deficit (which is of the order of the local 
macro-stress). Clearly the same result applies to those parts of the boundary that are 
parallel to the xj-axis. 

Next consider the case when the boundary may intersect the cells "diagonally". For a 
study of this case consider a simply connected domain c in B and let 5c be its boundary. 
Now consider the integral 

J (ajj(x) - Ejj(x))nj(x)ds 
5c 

where nj(x) denotes the outward unit normal to the boundary. 

It will prove useful in the sequel to observe from (2.21) and from the fact that Ijj(x) 
satisfies equilibrium that 

j (<7jj(x) - Ejj(x))n j(x)ds = j q(x)dA (2.30) 

5c c 

Choosing c as a square of size d and with sides parallel to the coordinate axes and by 
a repeated application of (2.28) using the fact that Hjjki(x) is d-periodic, it is found that 

J (djj(x) - Ey(x))nj(x)ds = plEjpJ 0(d 3 ) 

5c 

or equivalently (2.31) 
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j r,(x)dA = M-lEy uJ 0(d 3 ) 
c 

where the constant p. has the dimension of stress and c denotes a cell size and oriented 
domain and the bars indicate a suitable norm. Hence the residual is of order d^ but is also 
dependent on the second derivative of the strain field Ejj. Hence the errors are small 
provided the strain Ejj do not have larger gradients. 
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3 . APPT TCATTON TO PLATE WITH ORIENTED CRACKS 

The study by Litewka and Sawczuk [3] of the properties of plates perforated by 
periodically disposed slits is used as an example to study the effectiveness of the 
homogenization procedure. The example was selected because experimental data was 
available against which to test the predictions of the procedure. Furthermore, the shape and 
distribution of cracks is such that severe demands are made on the computations since the 
stress gradients within the unit elements are large in the vicinity of the slots. 

A representative drawing of the perforated plate is shown in Fig. 3.1a and the 
dimensions of the cracks are given in Fig. 3.1b. The plates were made from an aluminum 
alloy whose stress/strain relationship is given in Fig. 3.2. The macro stress, macro strain 
relationships for the perforated plate for different values of the crack alignment are also 
shown in Fig. 3.2. The failure modes are reported by Litwka and Morzynska [4] and 
involve various forms of strain locations as illustrated in Fig. 3.3. 


3.1 Trading and Boundary Conditions for Element 

Following the procedures described in Section 2, the periodic element is of size D = 
10 mm and the loading direction in relation to the direction of elongation of the cracks is 
defined by the angle a (Fig. 3.1b) The macro quantities are the stresses E x , I y , I xy and 
strains E x , E y , E xy while the micro quantities are represented by a x , a y , t xy and £ x , £ y , 
y xy repsectively. The homogenization procedure was completed using finite element 
representation of quarter of the element as shown in Fig. 3.5. All the analysis was 
performed using the ABAQUS [5] system using the 8 node, isoparametric plane stress 
element CPS8. Incremental elastic/plastic calculations were performed using the Mises 
Yield Criterion in conjunction with an isotropic incremental hardening law for which 


eR = - 
y 7 


f’(q) 

. a J 


sjjda 


where 


a 2 



(3.1) 
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and the uniaxial stress strain law has the form 


eP = f (o) 
df 

so that f'(a) = — 
da 

The function f(a) was obtained from the stress/strain relationship for the aluminum 
alloy given in Fig. 3.2. In this figure it should be noted that the stress reaches a limiting 
value of 155 MPa. 

The general case, subjected to macro stress £ x > ^y» ^xy> was analyzed using the 
element shown in Fig. 3.4. This covers the whole of the cyclically repeating cell. When 
the macro stress state is symmetric or anti- symmetric just a quarter of the repeating cell 
need to be analyzed (Fig. 3.5). 

Following the homogenization procedure described in Section 2 the displacement 
field has the form 


U(x,y) = U(x,y) + E x ■ x + E xy • y 
V(x,y) = V(x,y) + E y -y + E X y-x 

where U, V is the periodic component. Taking nodes A and D with the same value of x = a 
and nodes B and C with y = b. The boundary conditions imposed are; 


Ud * Ua “ Exy • D 
V D -V A = E yy D 

Uc ■ Ub = E X x ' E> 
v c - V B = E xy D 

The values of Exx, E yy and E xy are the loading inputs into the problem. 


(3.2) 


The general case 
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was analyzed under two types of loading: 


(i) strain controlled loading, where loading was applied by proportionately increasing the 
macroscopic strains E xx , E yy , E xy . 


(ii) stress controlled loading, when the loading was applied by proportionately increasing 
the resultant tractions F x , F y and F xy that are related to the macro stress 1^, ly and 
Zxy respectively. 


3.2 The Elastic Calculation 

The macro stress/strain relations in plane stress are written in the form: 


P _ Sxx _ V *y y 

xx m xx M yy yy 


y v 

p = 


J yy 


M yy M xx 


J yy 


y = 2E =^y 
i xy ^^xy 


G 


xy 


(3.3) 


From finite element calculation in the elastic range the following values for the 
homogenized elastic constants were determined: 

Elastic moduli M xx = 45700 MPa 
M yy = 64900 MPa 
G xy = 18700 MPa 

Poisson relation: v xy = 0.339 

v yx = 0.239 
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which satisfy the reciprocal relation 


= 

M y M x 


In matrix form the stress/strain relation in plane stress has the form 


^xx 


"21.9 -5.2 

0 " 


^XX 

Eyy 

= 10“ 6 

-5.2 15.4 

0 


lyy 

Yxy. 


0 0 

53.5_ 


^xy_ 



(3.4) 


For comparison with the constants given in Eq. (3.3) the value of Young's Modulus and 
the aluminum alloy of the plate is 69.5 MPa. It can be seen that the presence of the slots 
substantially reduce the stiffness of the system for shear and transverse loading. 


3.3 Elastic/Plastic Calculat ions 

The general case was analyzed under two types of loading. 

(i) Strain controlled loading, where loading was applied by proportionally 
increasing the strains E x , Ey and E X y. This form of loading resulted in stable 
calculations and is especially useful for determining the Limit Load Surface 
which is more fully described later. 

(ii) Stress controlled loading, where the loading was applied by proportionally 
increasing the stress £ x , I y and L xy . This is the form of loading reported in 
the tests of Litewka and Sawczuk. Convergence was rather poor when the 
loads approached the limit load value. 

Elastic/plastic calculations were performed for a uniaxial loading condition whose 
loading direction was defined by the angle a. This is the loading condition used by 
Litewke and Sawczuk. The calculations were repeated for different values of a which 
ranged between 0° and 90°. The computed macroscopic stress/strain results are shown in 
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Fig. 3.6 from which it may be seen that the experimental and calculated values are generally 
in good agreement. The predicted macroscopic stress/strain loading curves all show an 
initial elastic response, followed by an elastic/plastic region and a limit load condition for 
which the value of the load remains constant and further displacement is accommodated by 
plastic deformations. Since no failure criterion has been included in these calculations no 
failure is predicted from the calculations. In order to complete the relevant calculations a 
failure criterion must be included. In the absence of an obvious failure criterion this is not 
attempted at present, but shall be included later in an approximate analysis. 

The limit condition is a surface which is a function of Exy, Eyy and £xy The shape 
of this surface was obtained by performing elastic/plastic calculations for which the ratios 
of the strains remained constant. The calculations were repeated for various strain loading 
paths and the shape of the surfaces at the intersection with the section T*y = 0, section Z x 
= 0 and section X y = 0 are plotted in Figs. 3.7, 3.8 and 3.9 as are the strain rate vectors 
which satisfy the normality condition. The stress paths for different strain paths are shown 
in Fig. 3.7a and when the limit condition is reached the strain incremental vectors are 
normal to the surface. 

An isometric plot of the 3D yield surface is given in Fig. 3.10. Only the S xy > 0 
region is shown since the plot will be symmetrical with respect to the ^xx'^yy P^ ane - 

It is interesting to observe the shape of the Limit Surface for the intersection with Txy 
= 0. For this biaxial stress state the Limit Surface shown in Fig. 3.7 has the curious 
property that its shape is reminiscent of the Tresca yield condition in the first and third 
quadrants while the second and fourth quadrants indicate a condition more closely 
resembling the Mises condition. This topic is discussed more fully later. 

4. A STUDY OF EDGE EFFECTS 

In Section 2.3, it was illustrated how one of the errors introduced in the 
homogenization procedure is the imbalance of the applied macro stress and the 
calculated micro-stresses 6jj along the component boundary. For elastic conditions it was 
determined that the total imbalance force over the length of the boundary of length d is two 
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orders smaller than the value of the macro stresses. 

In order to study the influence of edge effects, an infinitely long strip in the x 
direction was analyzed. When the plate is subjected to axial stress L xx (Fig. 3.11) the 
sections such as AA and BB remain straight during deformation and consequently, only the 
shaded portion shown in Fig. 3.7 needs to be analyzed. 

The elastic solution for the case when = 0.02% and the elasto plastic solution for 
E xx = 1.5% have been determined and the variation of stress <7yy along the centerline is 
plotted as a function of the distance from the slot in Fig. 3.12. 

If the slotted system were analyzed at the macro-level then one of the boundary 
conditions would be Zyy = 0. In fact while Exx = 0 is satisfied on average the distribution 
of Cyy at the micro-level is the same as that shown in Fig. 3.12. It can be observed that the 
error in the edge stresses is given by Eq. (2.2) and when applied to the current problem 
gives 


pyy - Cyy * Eyy 


It is noted that 


D/2 

J Pyy dx = 0 
0 


D/2 

because j ayydx = 0 (Fig. 3.12a). 

0 

Hence the stresses obtained from the homogenized solution are in equilibrium with the 
applied macro-stress, provided equilibrium over the cell length is considered. However 
equilibrium is not satisfied at every point of the boundary. From St. Yenant principle it is 
to be expected that the homogenized solution is valid at a small distance from the edge. 
Reference to the solution of the strip problem indicate that the stresses and strains in the 
elements coincide everywhere with the homogenized solution with the exception of the 
edge elements. This observation is illustrated by Table 1 which gives the stress Cyy at 
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different locations in the strip. It can be seen that the stresses arc equal at the same location 
within the elements except for the edge elements, from which it can be deduced that the 
deviation from the homogenized solution is only one layer thick. This result is applicable 
for a continuous fiber reinforced composite loaded in the transverse direction but is not for 
loading in the longitudinal direction. 


Table 1 


Seven Cell Strip - a = 90 
Cyy Along Const Sections 
Elastic Solution - for E xx = 0.02% 


x 

y = 0.225 

10.225 

20.225 

30.225 

0.225 

-0.13 

-2.32 

-2.32 

-2.32 

1.000 

0.02 

-2.15 

-2.15 

-2.15 

1.775 

-0.09 

-1.96 

-1.95 

-1.95 

2.225 

-0.17 

-1.45 

-1.44 

-1.44 

3.000 

-0.42 

-0.25 

-0.24 

-0.24 

3.775 

0.30 

2.19 

2.19 

2.19 

4.056 

0.45 

3.28 

3.28 

3.28 

4.250 

0.31 

3.65 

3.64 

3.64 

4.444 

0.15 

3.92 

3.91 

3.91 

4.556 

0.42 

4.35 

4.34 

4.34 

4.750 

0.60 

4.66 

4.65 

4.65 

4.944 

0.39 

4.64 

4.63 

4.63 


The computations were extended to a loading level sufficiently high to cause 
extensive plastic deformations. The stress distribution Oyy along the centerline are shown 
in Fig. 12b which are in equilibrium with the macrostress lyy = 0. From the finite element 
calculations the observation can again be made that the stress distributions are the same for 
all elements except for those in the edge elements. This observation can be deduced from 
the results given in Table 2 which show the variation in the stress field Oyy along the x 
direction for the same section in adjacent elements. 
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Table 2 


Seven Cell Strip - a = 90 
Oyy Along y=const Sections 
Plastic Solution - E xx = 1.5% 


y = 0.225 


0.225 

1.000 

1.775 
2.225 
3.000 

3.775 
4.056 
4.250 
4.444 
4.556 
4.750 
4.944 


8.11 

4.36 

-12.30 

7.96 

15.10 

-6.55 

- 0.68 

-9.39 

-11.30 

-14.70 

- 2.00 

3.15 


10.225 


-66.70 

-52.90 

-42.70 

-29.30 

8.07 

39.00 

64.20 

79.60 

87.90 

94.50 

94.50 

94.70 


20.225 


-65.90 

-52.40 

-42.80 

-29.70 

7.21 

39.00 

64.40 

80.00 

88.40 
94.70 
94.50 
94.70 


30.225 


-65.90 

-52.40 

-42.80 

-29.70 

7.21 

39.00 

64.40 

80.00 

88.40 
94.70 
94.50 
94.70 


5. COMPUTATIONAL RESULTS AND COMPARISON WITH EXPERIMENT 

The predictions of the homogenization calculations are compared with the 
experimental investigation performed by Litewka and Sawczuk [1] in Table 3 with elastic 
properties defined by the relation (3.3) 


Table 3 



M x 

My 

G xy 

Uxy 

u yx 

Experiments 

41870 

60780 

19520 

0.294 

0.206 

Calculated 

45700 

64900 

18700 

0.339 

0.239 


The relations between the predictions and measurement are generally within 10%. 
The reciprocal relation 
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^xy _ ^yx 

M y M x 

is satisfied for the calculated constants and is closely satisfied by the experimental data 
which give the result 

o YV m 

— = 0.983 

Uyx My 

The macroscopic stress/strain relationships when the material is subject to uniaxial 
loading is given in Fig. 3.6 for different slot angles. The predictions compare with the 
experimental data with almost perfect agreement for all angle orientations except for a = 
7 t /6 when the difference in stress increases with strain from 0% to 10%. This may be due 
to early localization of the deformation for this orientation. The calculations are unable to 
predict the strains at which failure occurs. The large displacement analysis required to 
determine the failure strains has not been performed. Instead a simplified analysis of this 
problem has been completed and is described in the Appendix C 

The results of the approximate finite strain calculations are given in Figs. 3.13 and 
3.14 in which the stress and strain at failure are plotted for different slot angles a and strain 
hardening exponents. The results are normalized with respect to the failure stress and 
strain of an unperforated plate. The experimental failure stress values are also plotted in 
Fig. 3.13 and indicate generally good agreement with the predictions of the approximate 
analysis for low hardening exponents. The predictions of the failure strains in Fig. 3.14 
indicate that when a is greater than 20 ° there is a sudden decrease in the strain to failure. 
This prediction is in keeping with the experimental observations (Fig. 3.6). 

6 . DF.SCRTPTION OF LIMIT SURFACE AND CONSTITUTIVE EQUATION 

The isometric projection of the limit surface is shown in Fig. 3.10. The lower bound 
calculations performed in the Appendix define a limit surface defined by intersecting 
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surfaces with the form 


a ijk ilijl k l = °n (i,j,k,l = x,y) 

where ajj^l are defined in Eqn (B2) and a y is the material limit stress. There are eleven 
such relationships which indicate the complexity of the shape of the limit surface. In an 
attempt to simplify the problem an approximate analysis is developed in the Appendix B 
from which a simplified form of limit surface is given by three relationships instead of 

eleven 


n(z 


xx> Zyy-Zxy 


)■ 


a 


2 

n 


where 



The homogenization calculations have been used to develop the limit surface, but the 
means to formulate a closed-form mathematical description of the surface is less apparent. 
It is easy to fit Hill's generalized yield surface but this description does not pick up the 
shape shown in Fig. 3.7 for example. The procedure adopted here was to develop 
expressions for the surface based on approximate models and to select the constants of the 
approximate analysis to fit the results obtained from the homogenization procedures. 

It is also difficult to develop constitutive equations for the plate which are valid 
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current procedure is that computational experiments can be performed to determine the 
properties of the system. One such set of experiments consists of applying to the plate 
increments of strain such that the ratio 


AE 

AE 


M = 

yy 


k 


where k is a constant. The history of stress for some such strain paths is shown in Fig. 
(3.7). The initial response is elastic when the stress increments are dictated by the elastic 
constitutive equation of Eqn. (3.4). However as the applied stress approaches the limit 
surface the trajectory of the stress path is influenced increasingly by plasticity. The final 
location of the stress point on the limit surface is dictated by the normality rule when the 
vector of strain increments is normal to the limit surface. For some strain paths the stress 
state can change rapidly as the response of the material changes from purely elastic to 
purely plastic. This is particularly true of the paths 3 and 4 shown in Fig. 3.7. The 
simplest means of describing such behavior is to assume the component behaves in an 
elastic-perfectly plastic manner. With this assumption the stress paths would initially 
follow elastic response until the limit surface is reached when the normality rule would 
dictate the location of the stress state. Such an assumption gives stress histories which are 
close to those observed by calculation as indicated by the simple paths shown in Fig. 3.7. 
Then the simple constitutive relations are given by the sum of elastic deformation and 
plastic deformation 


^ij — CjjklSkl +Ejj 


where the matrix Qjki is given in Table 3 and the plastic deformation rate satisfies 


J jj 
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where C> has the form of the limit surface defined by Eqn. (6.1). 


CONCLUDING REMARKS 

It has been demonstrated that the procedures of homogenization can be used 
efficiently to determine the macroscopic properties of composite systems from the 
constituent properties. To investigate the efficiency of the procedure, a difficult problem 
was selected which involves very high stress gradients more severe than would normally 
be found in practice. The macroscopic response is highly anisotropic which is highly 
representative of the anisotropic effects of damage. Furthermore, both elastic and plastic 
properties have been investigated so that the application to practical systems with inelastic 
deformation should be able to proceed without difficulty. The internal state of composites 
is extremely complicated which underlines the need to formulate macroscopic constitutive 
equations with a limited number of state variables which represent the internal state at the 
micro level. 

As a next step it is intended to apply the method to determine the macroscopic 
properties of metal-matrix composites when subjected to the type of loadings to. be 
considered in the design process. The computational predictions are to be compared to the 
results of an experimental program. From this investigation constitutive equations shall be 
determined and couched in a form suitable for finite element calculations. An effort will be 
directed towards selecting the microstructure to optimize the performance during operation 
and the manufacturing process. 
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Some approximate analyses are performed to serve as a check on the validity of the 
homogenization procedures and to provide insight into the physical behavior of the plate. 

Approximate analyses for the elastic limit load and failure strains are now performed. 
These approximate calculations serve as check on the effectiveness of the homogeneization 
calculations and can also help to provide physical insight into the component behavior. 

6.1 Elastic Constants 

The elastic constants determined from the homogenization calculation are given in 
Table 3 when the elastic constants of the homogeneous material of the plate are 

E = 69550 MPa and o = 0.337 

Bounding estimates for E x and E y can be obtained using guessed equilibrium stress and 
compatable strain fields. 

To obtain an upper bound on M x a compatible strain field is assumed to be 
£ X = Ey = Ex 

where E x is the macrostrain. Energy balance then gives 

|M x e5 = i[|-e x 2 (0.95V) 

From these calculations the following results are obtained 
M x = 0.95E and o X y = 0.337 

To obtain a lower bound on M x the following stress field is assumed in the ligaments 
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between the slots 


<?x ~ » 


Gy =0 


Again using complementary energy balance 


1 Zt 1 


2 M x v 


= irj_(25: x )2ll 

vL2E x 


from which 


M x = 0.5E 

Similar procedures applied in the y direction give for the assumed strain calculations yields 
E y = 0.95E 

while the stress field calculation gives 


E y = 0.9E 

Calculations have also been performed by assuming the slots do not interact and the energy 
icnrease in the body is given by the energy change for a slot in an infinite body. 

For calculations in the x direction and assuming that the slots can be represented by 
slits the elastic internal energy per unit area of the plate is found to be 


z\ jta 2 Z x 

— *- + -n 

2E 2E 2E 
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where n is the number of cracks per unit area and the cracks are assumed to be non- 
interacting. 

In the present problem a = 2.5mm and n = 1/100 so that 


E x = 


l + 7ta n 


0.72E 


The summary of results of experimental and computed and approximate results for elastic 
constants. 


Lower 



Non-Interacting 

Upper 

Bound 

Experiment 

Horn 

Cracks 

Bound 

x direction 0.5E 

0.602E 

0.66E 

0.72E 

0.95E 

y direction 0.9E 

0.874E 

0.933E 

1.00E 

0.95E 


The predictions of the homogenization procedure agree to within 10% of the 
experimental observations. The non-interacting results are the most consistent of the 
approximate procedures which use extremely simple calculations. 
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APPENDIX B - T-TMTT LOAD SURFACE 

The isometric projection of the limit surface obtained from the homogenization 
calculations is shown in Fig. 3.10. In spite of the fact that the Mises yield condition is 
used it can be seen that the yield surface for the plate demonstrates features of both Tresca 
and Mises type criteria. In order to study the shape of the yield surface more fully 
approximate solutions are now developed. 

T x~>wer Bound Stress Analysis 

An equilibrium stress field for the plate is illustrated in Fig. A.l. It can be deduced 
that there are nine microstress fields o xx , ^yy> CT xy gi ven 1° terms of the macro stresses 
I xx> Eyy and E xy applied to the perimeter of the plate element. In terms of the 
macrostress E^ Ey, 2^y the element perimeter stresses are given by 

E xx = 0.5a x Eyy = 0.9a y E xy = 0.4x (B.l) 

The value of o x , a y and x for the nine stress fields are given in Table A.l and when these 
are substituted in the yield condition ct xx + c yy - <J X x a yy + 3a xy = a 2 then the following 
yield criteria are established. 

Field 1) a y + 3x 2 < c 2 

n O v T 9 9 

Field 2) + -^- + x z < c£ 

£ 

Field 3) a x + a y - c x c y + (a x + a y )^ + x 2 < a 2 

Field 4) aj < a 2 

Field 5) a 2 <a 2 ( B - 2 ) 
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Field 6) 

_2 , £££ , _2< _2 
a x + ^ +x sa 0 



Field 7) 

a x + a y “ a x CT y + 3x2 - a o 



Field 8) 

o x + 3x 2 < a 2 



Field 9) 

a x + a 2 -a x a y <a 2 



Table A 1 




Field 

a xx 

a yy 

a xy 

1 

0 

a y 

X 

2 

t/2 

Cy + X/2 

x/2 

3 

<t x + t/2 

Gy +X/2 

x/2 

4 

<*x 

0 

0 

5 

0 

ay 

0 

6 

a x + x/2 

x/2 

x/2 

7 

°x 

ay 

X 

8 

^x 

0 

X 

9 

o x 

a y 

0 


The inequalities given in Eqns. (B.2) define lower bounds of the limit surface. 
Sections of the limit surface and the results of the lower bound calculations are shown in 
Figs. A.2-A.3 and A.4. Also shown in the figures are the yield conditions for the 
undamaged system. As is to be expected the equilibrium stress fields give a lower bound 
on the yield surface. Generally, the bounds are good except in the Eyy, I xy intersection 
where the difference between the lower bound and the homogenization results can be 
considerable. 
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Approximate Limit Surface An alysis 

Another analysis is performed using engineering approximations but which give 
results agreeing more closely with the numerical results. 

Referring to Fig. A.5 the plate element has been divided into three sections I, II and 
in which represents the critical regions where unconstrained plastic deformation can occur. 
The element I is subjected to the stress state 


CT n - 


T nt ~ 2Lxy 


With the assumption that the strain in the t direction is e t = 0 gives the result 


a t = — = Z x 

i 2 x 


Using this stress state in the limit condition 




where a u is the ultimate stress of the material gives the condition 


/V \ 2 fr. \ 2 




+ 12 


■'xy 

V°uy 


= 1 


(B.l) 


The element H is subjected to the average stress state 


v f 10 

° n = Ht 


-T f 10 

°nt “ ^xy^ g 


With the assumption that the stress Of = 0 the limit load criterion becomes 
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(B.2) 


1 = 




With a stress field on element m given by 


<y n — 2L X 



(B.3) 


Then the limit condition for this element is 


or 


cJ = aS + a?-a n a t + 3a 


2 

tn 



K°n) 



(z. 






20 S x S y 

s h J — 

9 u; 


fz S 2 

^xy 
^ a n ; 


(B.4) 


The estimated yield conditions (B.1)-(B.3) are plotted in Figs. A. 2, A. 3 and A. 4 
which also show the calculated limit conditions. The conditions (B.l) and (B.2) compare 
well with the calculated values. However the predictions of (B.3) are less good. 

To improve condition (B.4) it is modified by adjusting the constants in (B.3) so that it 
fits closely the relationship determined from the homogenization calculations. The 
relationship (B.4) is modified to give 


^ v nocrp ^ 




V 81 ) 


v CT uy 


- 1.2 


( . 19 


v°uy 


a J + 4 


_ 


= 1 


V a u 


(B.5) 


A simplified descripton of the limit surface is then given by the surfaces 


& 


^2 


+ 12 


uj 


‘'xy 

V a uy 


= 1 
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100 

81 






( 

2' 

+ 3 

^xy 






= 1 


8 

3 



+ 


lOO ^j 

81 J 




V°u7 


- 1.2 




v°uA°u; 




K<*u) 


(B.6) 
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APPENDIX C: APPROXIMATE A NALYSIS OF FAILURE STRAIN S 

An approximate finite strain analysis is now performed when the perforated plate is 
subject to uniaxial stress for different slot orientations. The failure analyses are applied to 
the two characteristic ligaments which experiments indicate are the sources of failure. 
Referring to Fig. A.6 the ligment I is the ligament in line with adjacent slots while ligament 
II is the ligament between laterally adjacent slots. 

The constitutive law used is 

a e = Ke" 

where a, is the effective stress, e e the effective strain and K and n are material constants. 

V 

The strain components are defined by the relationship 



where Sjj is the stress deviator. 

The classical stability analysis for an undamaged homogeneous material subjected to a 
uniaxial stress state predicts the failure strain at the load maximum to be 


e° = n 

and the nominal stress at load maximum is 
a° = Kn n exp(-n) 

Reference to Fig. 3.6 indicates a failure strain of the undamaged material to be 1.5% so that 
n for the material in question is 
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n = 0.015 


Analysis of Ligament I 

The plate is subjected to the macro stress a shown in Fig. A. 6 with stresses in 
ligament I 

2A n D . 2 2A 0 D 

a n = — 2 o/ sin a x nt = — — sinacosa 

n A ( D-£y nt A 1 D-/y 

where A 0 is the original cross sectional area and A is the current cross sectional area of the 
ligament. 

Making the assumption that the strain in the t direction is zero i.e. £[ = 0 gives 



With the given approximate stress state the effective stress a e is 


a e = V3sina 


1-— sin 2 a 
4 


nl/2 


a, 


D 


A *D-/y 


Assuming a deformation theory for the plastic strains the ratio for e n /e nt can be determined 
in terms of stress deviators so that 


e n _ s n _ 1 2o n - o t _ 1 sing 


^nt s nt ^ ^ 


nt 


2 cos a 


and the effective strain is then 
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and 


2 / 2 \ 1 / 2 

8 e = -^=^l + 4 c °t aj e n 


From the definition of finite strain 


— = exp(-e n ) 


so that the stress/normal ligament strain e n relationship becomes 


c U = 


n 


2/_l ^l + 4cot Z aj 




V3 


sin a 


- — sin 2 a 

4 J 


1/2 Ke n exp( e n )- D 


From this expression the strain at the load maximum following the classical analysis is 


e n = n. 


The strain in the longitudinal direction in the band is 


ej = e n sin 2 a + 2e nt sinacosa 


and recalling that 


cos a 


e nt = 2-r— ' e n 


sin a 


yields the strain in the longitudinal direction as 


4= 


sin a + 4cos‘ 


! a]e n 
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Referring to Fig. A.6 the average strain is 


(^/sintt^ ^r 2 ; | 

* (D/sina) D !■ J n 


At load maximum e n = n so that 


= — [sin 2 a + 4cos 2 a] 

e? d l j 


is the average failure strain. 


Analysis of Ligament II 

A similar analysis for ligament follows along similar lines except that the stress a t is 
taken to be zero. Then 

A 0 2 D 

a n = — ^ Ot cos a 

A D-* x 

A 0 D 

x nt = — ^a^sinacosa 

A D -£ x 

a t = 0 

Hence the effective stress a e is 

a e = cosa [3 - 2cos 2 a]^ 2 — 5 — 

e A D-* x 

and the corresponding strains are 
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e n _ 2cosa 
e nt 3sina 


e? = [l + 3tan 2 a]£ 2 


The stress Of — £ n relationship is then given by 


[l + 3tan 2 a] n/2 ^ n , 

Gf = — 5 — 7172 ' K-En ex Pt e n/ ^ 

* cosa[3-2cos 2 a] 1/2 D 


In the band the strain in the longitudinal direction is 


E^ 1 = [sin 2 a + 4cos^a]£ n 


and the average strain in the longitudinal direction is 


«y/sina) e£ . ty [si „2 a + 4cos 2 a]e 
* (D/sina) D 


At the load maximum £ n = n which gives 


E, = — [sin 2 a + 4cos 2 a]n 
4 D 


and 


— = ^[sin 2 cx + 4cos 2 a] 

£? D L 
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Fig. 2.1 A structure made of a periodic composite material. 
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Fig. 2.2a A nine-cell sub-domain within the body shown in Fig. 2.1 
Pj, P 2 . . .: homologous points. ABCD: unit cell. 
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Fig. 2.2b Deformed shape of nine-cell sub-domain; 

» 1 , R radius of circle passing through three consecutive homologous points. 
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Fig. 2.2c Macroscopically homogenous deformation of nine-cell 
sub-domain is a good approximation to the deformed 
shape of Fig. 2.2b. 
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The program of homogenization fails locally if the 
deformed shape sub-domain of the homogenized body 
is such that R is comparable to d. Such a situation 
can develop near a crack tip or when localized defor- 
mation occurs* 





Fig. 3.1 


a) Specimen with simulated damage. 

b) Dimensions of cracks. 
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Fig. 3.2 Macroscopic stress-strain curves for specimens with 
different crack orientations. 
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Fig. 3.3 Strain localization at failure. 
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^8* 3.5 Finite element representation used in homogenization process. 
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Fig. 3.6 Comparison of computed and experimental 

results. 
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Fig. 3.8 Limit Surface in I I^y plane. 
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Fig. 3.9 Limit Surface in Ey y , plane. 
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Fig. 3.10 Isometric view of limit surface. 
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Fig. A. 1 Equilibrium stress fields in element. 
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Fig. A.2 Comparison of Calculation and Approximate Limit Surface 

in I xx , I yy plane. 
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Fig. A. 3 Comparison of Calculation and Approximate Limit Surface 

in Z xx , plane. 
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Fig. A.4 Comparison of Calculation and Approximate Limit Surface 

in Zyy, Zxy plane. 
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Fig. A. 6 Model for simplified analysis of stress and strain at load maximum. 
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